function I = hwozd_Image(Im)
%UNTITLED1 Summary of this function goes here
%  Detailed explanation goes here
I.image = im2double(Im);
%I.x0 = [];
%I.y0 = [];
%I.x = [];
%I.y = [];
%I.m00 = -1;

    d = size(I.image);
    xstep		= 2/(d(1)-1);
    ystep		= 2/(d(2)-1);
    [x,y]		= meshgrid(-1:ystep:1,-1:xstep:1);

    %----------	compute unit disc mask
    circle1		= x.^2 + y.^2;
    inside		= find(circle1<=1);
    mask		= zeros(d);
    mask(inside)	= ones(size(inside));
    ind = find(mask == 0);
    x(ind) = 0;
    y(ind) = 0;
    I.x = x;
    I.y = y;
    %I.x0int =
    
    %I.m00 = hwozd_m(I,0,0);
    xs = ones(size(I.x));           %(I.x) .^0;
    ys = ones(size(I.y));           %(I.y) .^0;
    I.m00 = sum(sum( (xs .* ys) .* I.image));
    
    %x0int
    p = 1;
    q = 0;
    xs = repmat((1:d(1))',1,d(2));
    ys = repmat((1:d(2)),d(1),1);
    xs = xs .^p;
    ys = ys .^q;
    moment = sum(sum( (xs .* ys) .* double(I.image)));
    I.x0int = round(moment / I.m00);
    
    %y0int
    p = 0;
    q = 1;
    xs = repmat((1:d(1))',1,d(2));
    ys = repmat((1:d(2)),d(1),1);
    xs = xs .^p;
    ys = ys .^q;
    moment = sum(sum( (xs .* ys) .* double(I.image)));
    I.y0int = round(moment / I.m00);
    
    
    %x0 = hwozd_m(I,1,0) / I.m00;
    xs = (I.x);                         %(I.x).^1
    ys = ones(size(I.y));               %(I.y) .^0;
    x0 = sum(sum( (xs .* ys) .* I.image)) / I.m00;
    
    %y0 = hwozd_m(I,0,1) / I.m00;
    xs = ones(size(I.x));           %(I.x).^0;
    ys = (I.y);                     %(I.y).^1
    y0 = sum(sum( (xs .* ys) .* I.image)) / I.m00;
    
    x0_ = zeros(d) + x0;
    y0_ = zeros(d) + y0;
    
    x0_(ind) = 0;
    y0_(ind) = 0;
    
    I.x0 = x0_;
    I.y0 = y0_;
    
    
I = class(I,'hwozd_Image');